rm(list = ls())
#working_folder<- "/Users/bgpopescu/Library/CloudStorage/Dropbox/covid_institutions/"
working_folder<-"//Mac/Dropbox-1/covid_paper_replication/"
setwd(working_folder)

#This is to obtain:
#-figure_a10a
#-figure_a10b
#-figure_a10c
#-figure_a10d
#-figure_a10e
#-figure_a10f
#-figure_a10g
#-figure_a10h
#-figure_a10i
#-figure_a10j
#-figure_a10k


library("readxl")
library("reshape")
library("stringi")
library("stringr")
library("dplyr")
library("tidyr")
library("stargazer")
library("corrplot")
library("Hmisc")
library("glue")
library("ggplot2")
library('gridExtra')
library('grid')
library('lfe')
library('broom')
library("plm")
library("multiwayvcov")
library("lmtest")
library("psych")
library("ggrepel")
library("lemon")
library("sandwich")
library("spdep")
library("xtable")
library('gsynth')
library("pracma")
library("conleyreg")
library("vars")
library("rlist")
library("pgirmess")
library("geojsonsf")



#Step1. Turning some variables to numeric
data_cov_mob = geojson_sf("./data/data_cov_mob.geojson")
data_cov_mob$pct_manufacturing<-as.numeric(data_cov_mob$pct_manufacturing)

data_cov_mob$kul<-as.numeric(data_cov_mob$kul)
data_cov_mob$nat<-as.numeric(data_cov_mob$nat)
data_cov_mob$frei<-as.numeric(data_cov_mob$frei)
data_cov_mob$soz<-as.numeric(data_cov_mob$soz)
data_cov_mob$pol<-as.numeric(data_cov_mob$pol)
data_cov_mob$inter<-as.numeric(data_cov_mob$inter)

data_nocovid_nomob<-subset(data_cov_mob, week=="00")

data_nocovid_nomob$kul<-as.numeric(data_nocovid_nomob$kul)
data_nocovid_nomob$nat<-as.numeric(data_nocovid_nomob$nat)
data_nocovid_nomob$frei<-as.numeric(data_nocovid_nomob$frei)
data_nocovid_nomob$soz<-as.numeric(data_nocovid_nomob$soz)
data_nocovid_nomob$pol<-as.numeric(data_nocovid_nomob$pol)
data_nocovid_nomob$inter<-as.numeric(data_nocovid_nomob$inter)

data_cov_mob$gspo<-as.numeric(data_cov_mob$gspo)
data_cov_mob$gkul<-as.numeric(data_cov_mob$gkul)
data_cov_mob$gnat<-as.numeric(data_cov_mob$gnat)

data_cov_mob$gsoz<-as.numeric(data_cov_mob$gsoz)
data_cov_mob$gpol<-as.numeric(data_cov_mob$gpol)
data_cov_mob$ginter<-as.numeric(data_cov_mob$ginter)

#Bridging Networks
data_cov_mob$bridging_total<-data_cov_mob$gspo+ data_cov_mob$gkul+ data_cov_mob$gnat
data_cov_mob$bridging_tot_pop<-data_cov_mob$bridging_total/data_cov_mob$pop_total*1000
mean_brid<-summary(data_cov_mob$bridging_tot_pop)[4]
data_cov_mob$bridging_tot_pop_dummy<-ifelse(data_cov_mob$bridging_tot_pop>mean_brid, 1, 0)

#Bonding Networks
data_cov_mob$bonding_total<-data_cov_mob$gsoz+ data_cov_mob$gpol+ data_cov_mob$ginter
data_cov_mob$bonding_tot_pop<-data_cov_mob$bonding_total/data_cov_mob$pop_total*1000
mean_bond<-summary(data_cov_mob$bonding_tot_pop)[4]
data_cov_mob$bonding_tot_pop_dummy<-ifelse(data_cov_mob$bonding_tot_pop>mean_bond, 1, 0)

#All Networks
data_cov_mob$gvereine<-as.numeric(data_cov_mob$gvereine)
data_cov_mob$vereine_tot_pop<-data_cov_mob$gvereine/data_cov_mob$pop_total*1000
mean_ver<-summary(data_cov_mob$vereine_tot_pop)[4]
data_cov_mob$vereine_tot_pop_dummy<-ifelse(data_cov_mob$vereine_tot_pop>mean_ver, 1, 0)

data_cov_mob_week<-subset(data_cov_mob, week=="01")
data_cov_mob_week<-subset(data_cov_mob_week, select = c("GEN", 
                                                        "vereine_tot_pop",
                                                        "bridging_tot_pop", "bonding_tot_pop"))

data<-data_cov_mob
data$week_new<-as.numeric(data$week)+1
data$week_new<-paste0("0", data$week_new)  
data$week_new2<-ifelse(nchar(data$week_new)==3, str_remove(data$week_new, "^0+"), data$week_new)
data$week<-data$week_new2
data<-subset(data, select = -c(week_new, week_new2))
data$post_treat<-ifelse(as.numeric(data$week)>10, 1,0)


data$post_treat_vereine_tot_pop_dummy<-data$post_treat*data$vereine_tot_pop_dummy
data$post_treat_bridging_tot_pop_dummy<-data$post_treat*data$bridging_tot_pop_dummy
data$post_treat_bonding_tot_pop_dummy<-data$post_treat*data$bonding_tot_pop_dummy
data$post_treat_pct_service_sec_narrow<-data$post_treat*data$pct_service_sec_narrow
data$post_treat_change_afd<-data$post_treat*data$pct_change_afd
data$post_treat_police_stations_area<-data$post_treat*data$police_stations_area

#############################
#Function for calculating SE#
#############################

# Calculate robust confidence intervals
se_robust_cluster <- function(x){
  res<-coeftest(x, vcov = cluster.vcov(x, data$kreis_name, force_posdef=T))[, 2]
  return(res)}

ci_robust_cluster <- function(x){
  res_ci<-coefci(x, vcov = cluster.vcov(x, data$kreis_name, force_posdef=T))
  return(res_ci)}


###########
#All clubs#
###########


#Create dummies

library('fastDummies')
# Create dummy variables:
dataf <- dummy_cols(data, select_columns = 'AGS')
dataw <- dummy_cols(dataf, select_columns = 'week')
datas <- dummy_cols(dataw, select_columns = 'SN_L')

library(dplyr)
new<-dataw %>% dplyr::select(starts_with("AGS_"))
new2<-dataw %>% dplyr::select(starts_with("week_"))
new3<-datas %>% dplyr::select(starts_with("SN_L"))

list_dummies<-names(new)
#I am removing the counties which get dropped in Stata

#Re-enable this later if necessary
list_dummies<-list_dummies[list_dummies!= "AGS_0"]
list_dummies<-list_dummies[list_dummies!= "AGS_5316"]
list_dummies<-list_dummies[list_dummies!= "AGS_8121"]
list_dummies<-list_dummies[list_dummies!= "AGS_8221"]
list_dummies<-list_dummies[list_dummies!= "AGS_8231"]
list_dummies<-list_dummies[list_dummies!= "AGS_8311"]
list_dummies<-list_dummies[list_dummies!= "AGS_8421"]
list_dummies<-list_dummies[list_dummies!= "AGS_9461"]
list_dummies<-list_dummies[list_dummies!= "AGS_13004"]
list_dummies<-list_dummies[list_dummies!= "AGS_16052"]
list_dummies<-list_dummies[list_dummies!= "AGS_16053"]
list_dummies<-list_dummies[list_dummies!= "AGS_16054"]
list_dummies<-list_dummies[list_dummies!= "AGS_16055"]
list_dummies<-list_dummies[list_dummies!= "AGS_16061"]
list_dummies<-list_dummies[list_dummies!= "AGS_16077"]

list_week_dummies<-names(new2)
list_week_dummies<-list_week_dummies[list_week_dummies!= "week_new"]
list_week_dummies<-list_week_dummies[list_week_dummies!= "week_new2"]
list_week_dummies<-list_week_dummies[list_week_dummies!= "week_numeric"]
#Re-enable this later if necessary
list_week_dummies2<-list_week_dummies[list_week_dummies!= "week_09"]
list_week_dummies2<-list_week_dummies2[list_week_dummies2!= "week_52"]


list_land_dummies<-names(new3)
#Re-enable this later if necessary
list_land_dummies_full<-list_land_dummies[list_land_dummies!= "SN_L"]
list_land_dummies<-list_land_dummies_full[list_land_dummies_full!= "SN_L_01"]


#Event study
#Bridging Interaction
new2<-dataw %>% dplyr::select(starts_with("week_"))
list_week_dummies<-names(new2)
list_week_dummies<-list_week_dummies[list_week_dummies!= "week_new"]
list_week_dummies<-list_week_dummies[list_week_dummies!= "week_new2"]
list_week_dummies<-list_week_dummies[list_week_dummies!= "week_numeric"]

#######################################
#Interaction terms for all indep. vars#
#######################################

#Creating the interaction terms for bridging
for (i in list_week_dummies){
  datas[[paste("bri_", i, sep="")]]<-datas[[i]]*datas$bridging_tot_pop_dummy
  datas[[paste("pct_turnout_", i, sep="")]]<-datas[[i]]*datas$pct_turnout
  datas[[paste("pop_total_", i, sep="")]]<-datas[[i]]*datas$pop_total
  datas[[paste("log_gdp_", i, sep="")]]<-datas[[i]]*datas$log_gdp_per_cap
  datas[[paste("pop_den_", i, sep="")]]<-datas[[i]]*datas$pop_den
  datas[[paste("pct_college_", i, sep="")]]<-datas[[i]]*datas$pct_college
  datas[[paste("pct_above_60_", i, sep="")]]<-datas[[i]]*datas$pct_pop_above_60
  datas[[paste("pct_sec_", i, sep="")]]<-datas[[i]]*datas$pct_service_sec_narrow
  datas[[paste("pct_manuf_", i, sep="")]]<-datas[[i]]*datas$pct_manufacturing
  datas[[paste("police_den_", i, sep="")]]<-datas[[i]]*datas$police_stations_area
  datas[[paste("east_ger_", i, sep="")]]<-datas[[i]]*datas$east_ger
}

#Step2 Selecting the right independent vars
inter_terms_bri<-datas %>% dplyr::select(starts_with("bri_week"))
inter_terms_bri<-names(inter_terms_bri)
inter_terms_bri<-inter_terms_bri[inter_terms_bri != "bri_week_10"]

inter_terms_pct_turnout<-datas %>% dplyr::select(starts_with("pct_turnout_week"))
inter_terms_pct_turnout<-names(inter_terms_pct_turnout)
inter_terms_pct_turnout<-inter_terms_pct_turnout[inter_terms_pct_turnout != "pct_turnout_week_10"]

inter_terms_pop_total<-datas %>% dplyr::select(starts_with("pop_total_week"))
inter_terms_pop_total<-names(inter_terms_pop_total)
inter_terms_pop_total<-inter_terms_pop_total[inter_terms_pop_total != "pop_total_week_10"]

inter_terms_log_gdp<-datas %>% dplyr::select(starts_with("log_gdp_week"))
inter_terms_log_gdp<-names(inter_terms_log_gdp)
inter_terms_log_gdp<-inter_terms_log_gdp[inter_terms_log_gdp != "log_gdp_week_10"]

inter_terms_pop_den<-datas %>% dplyr::select(starts_with("pop_den_week"))
inter_terms_pop_den<-names(inter_terms_pop_den)
inter_terms_pop_den<-inter_terms_pop_den[inter_terms_pop_den != "pop_den_week_10"]

inter_terms_pct_college<-datas %>% dplyr::select(starts_with("pct_college_week"))
inter_terms_pct_college<-names(inter_terms_pct_college)
inter_terms_pct_college<-inter_terms_pct_college[inter_terms_pct_college != "pct_college_week_10"]

inter_terms_pct_above_60<-datas %>% dplyr::select(starts_with("pct_above_60_week"))
inter_terms_pct_above_60<-names(inter_terms_pct_above_60)
inter_terms_pct_above_60<-inter_terms_pct_above_60[inter_terms_pct_above_60 != "pct_above_60_week_10"]

inter_terms_pct_sec<-datas %>% dplyr::select(starts_with("pct_sec_week"))
inter_terms_pct_sec<-names(inter_terms_pct_sec)
inter_terms_pct_sec<-inter_terms_pct_sec[inter_terms_pct_sec != "pct_sec_week_10"]

inter_terms_pct_manuf<-datas %>% dplyr::select(starts_with("pct_manuf"))
inter_terms_pct_manuf<-names(inter_terms_pct_manuf)
inter_terms_pct_manuf<-inter_terms_pct_manuf[inter_terms_pct_manuf != "pct_manuf_week_10"]

inter_terms_police_den<-datas %>% dplyr::select(starts_with("police_den"))
inter_terms_police_den<-names(inter_terms_police_den)
inter_terms_police_den<-inter_terms_police_den[inter_terms_police_den != "police_den_week_10"]

inter_terms_east_ger<-datas %>% dplyr::select(starts_with("east_ger"))
inter_terms_east_ger<-names(inter_terms_east_ger)
inter_terms_east_ger<-inter_terms_east_ger[inter_terms_east_ger != "east_ger_week_10"]


#Deleting variables which are collinear
list_land_dummies_full2<-list_land_dummies_full[list_land_dummies_full!= "SN_L_02"]
list_week_dummies2<-list_week_dummies[list_week_dummies!= "week_53"]

#Step3 Creating the list of variables
for_event_inter_vars_complete<-list.append(inter_terms_bri,
                                                inter_terms_pct_turnout,
                                                inter_terms_pop_total,
                                                inter_terms_log_gdp,
                                                inter_terms_pop_den,
                                                inter_terms_pct_college,
                                                inter_terms_pct_above_60,
                                                inter_terms_pct_sec,
                                                inter_terms_pct_manuf, 
                                                inter_terms_police_den,
                                                inter_terms_east_ger,
                                                list_land_dummies_full2, 
                                                list_week_dummies2)

non_int_vars<-c("lag_covid_pc",
  "pct_turnout",
  "log_pop_total",
  "log_gdp_per_cap",
  "pop_den",
  "pct_college",
  "pct_pop_above_60",
  "pct_service_sec_narrow",
  "pct_manufacturing",
  "east_ger")

for_event_inter_vars1<-list.append(inter_terms_bri,
                                  inter_terms_pct_turnout,
                                  non_int_vars,
                                  list_land_dummies_full2, 
                                  list_week_dummies2)

for_event_inter_vars2<-list.append(inter_terms_bri,
                                   inter_terms_pop_total,
                                   non_int_vars,
                                   list_land_dummies_full2, 
                                   list_week_dummies2)

for_event_inter_vars3<-list.append(inter_terms_bri,
                                   inter_terms_log_gdp,
                                   non_int_vars,
                                   list_land_dummies_full2, 
                                   list_week_dummies2)

for_event_inter_vars4<-list.append(inter_terms_bri,
                                   inter_terms_pop_den,
                                   non_int_vars,
                                   list_land_dummies_full2, 
                                   list_week_dummies2)

for_event_inter_vars5<-list.append(inter_terms_bri,
                                   inter_terms_pct_college,
                                   non_int_vars,
                                   list_land_dummies_full2, 
                                   list_week_dummies2)

for_event_inter_vars6<-list.append(inter_terms_bri,
                                   inter_terms_pct_above_60,
                                   non_int_vars,
                                   list_land_dummies_full2, 
                                   list_week_dummies2)

for_event_inter_vars7<-list.append(inter_terms_bri,
                                   inter_terms_pct_sec,
                                   non_int_vars,
                                   list_land_dummies_full2, 
                                   list_week_dummies2)

for_event_inter_vars8<-list.append(inter_terms_bri,
                                   inter_terms_pct_manuf,
                                   non_int_vars,
                                   list_land_dummies_full2, 
                                   list_week_dummies2)

for_event_inter_vars9<-list.append(inter_terms_bri,
                                   inter_terms_police_den,
                                   non_int_vars,
                                   list_land_dummies_full2, 
                                   list_week_dummies2)


for_event_inter_vars10<-list.append(inter_terms_bri,
                                   inter_terms_east_ger,
                                   non_int_vars,
                                   list_land_dummies_full2, 
                                   list_week_dummies2)



#Step4 Deleting missing observations
datas$time<-as.numeric(datas$week)
for_event_inter_vars_complete<-append(for_event_inter_vars_complete,
                              c("lag_covid_pc",
                                "pct_turnout",
                                "log_pop_total",
                                "log_gdp_per_cap",
                                "pop_den",
                                "pct_college",
                                "pct_pop_above_60",
                                "pct_service_sec_narrow",
                                "pct_manufacturing",
                                "east_ger"))

for_event_inter_vars11<-list.append(for_event_inter_vars_complete, 
                                    list_land_dummies_full2)





all_relev_vars<-subset(datas, select = list.append(for_event_inter_vars_complete, 
                                                     "bridging_tot_pop_dummy", 
                                                     "mean_mobil", "POINT_X", "POINT_Y",
                                         "time", "AGS"))
all_relev_vars<-all_relev_vars[complete.cases(all_relev_vars), ]

#Step5 Deciding on the temporal lag - value with the lowest AIC
x<-VARselect(all_relev_vars$mean_mobil)
min(x$criteria["AIC(n)",])
temp_lag<-names(x$criteria["AIC(n)",])[x$criteria["AIC(n)",]==min(x$criteria["AIC(n)",])]
temp_lag<-as.numeric(temp_lag)

#Step6 Merging with spatial dataframe
nuts<-subset(data_cov_mob, week=="00")
nuts2<-subset(nuts, select = c("AGS"))
nuts2$AGS<-as.numeric(as.character(nuts2$AGS))
all_relev_vars_oneweek<-subset(all_relev_vars, week_11==1)


#Step7 Merging with shape file
nuts_merged<-left_join(nuts2, all_relev_vars_oneweek, by = c("AGS" = "AGS"))
nuts_merged2<-subset(nuts_merged, !is.na(nuts_merged$bridging_tot_pop_dummy))

#Step8 Redefinining projection to calculate dist in meters
nuts_merged3 <- st_transform(nuts_merged2, crs = "+proj=utm +zone=1 +units=cm")
nuts_merged3$POINT_X_m<-st_coordinates(st_centroid(nuts_merged3))[,1]
nuts_merged3$POINT_Y_m<-st_coordinates(st_centroid(nuts_merged3))[,2]

#Step9: Calculating the correlogram
coo<-subset(nuts_merged3, select = c("POINT_X_m", "POINT_Y_m"))
coo<-st_drop_geometry(coo)

#Calculating the correlogram
pgi_cor <- correlog(coords=coo, z=nuts_merged3$bridging_tot_pop_dummy, method="Moran", 
                    nbclass=10)
pgi_cor_df<-data.frame(pgi_cor)
dist_cut_all_vars<-round(pgi_cor_df$dist.class[abs(pgi_cor_df$coef)== min(abs(pgi_cor_df$coef))]/1000,0)


#Step10: pct_service formula
time_trend_all_vars<-as.formula(paste("mean_mobil~",
                                 paste(for_event_inter_vars_complete, collapse= "+")))
dm_all_relev_vars <- dist_mat(all_relev_vars, lat = "POINT_X", lon = "POINT_Y", unit = "AGS", time = "time")


time_trend_all_vars1<-as.formula(paste("mean_mobil~",
                                      paste(for_event_inter_vars1, collapse= "+")))

time_trend_all_vars2<-as.formula(paste("mean_mobil~",
                                       paste(for_event_inter_vars2, collapse= "+")))
time_trend_all_vars3<-as.formula(paste("mean_mobil~",
                                       paste(for_event_inter_vars3, collapse= "+")))
time_trend_all_vars4<-as.formula(paste("mean_mobil~",
                                       paste(for_event_inter_vars4, collapse= "+")))

time_trend_all_vars5<-as.formula(paste("mean_mobil~",
                                       paste(for_event_inter_vars5, collapse= "+")))

time_trend_all_vars6<-as.formula(paste("mean_mobil~",
                                       paste(for_event_inter_vars6, collapse= "+")))

time_trend_all_vars7<-as.formula(paste("mean_mobil~",
                                       paste(for_event_inter_vars7, collapse= "+")))

time_trend_all_vars8<-as.formula(paste("mean_mobil~",
                                       paste(for_event_inter_vars8, collapse= "+")))

time_trend_all_vars9<-as.formula(paste("mean_mobil~",
                                       paste(for_event_inter_vars9, collapse= "+")))


time_trend_all_vars10<-as.formula(paste("mean_mobil~",
                                       paste(for_event_inter_vars10, collapse= "+")))

time_trend_all_vars11<-as.formula(paste("mean_mobil~",
                                        paste(for_event_inter_vars11, collapse= "+")))


x<-lm(time_trend_all_vars, data=all_relev_vars)
summary(x)

x1<-lm(time_trend_all_vars1, data=all_relev_vars)
summary(x1)

x2<-lm(time_trend_all_vars2, data=all_relev_vars)
summary(x2)

x3<-lm(time_trend_all_vars3, data=all_relev_vars)
summary(x3)

x4<-lm(time_trend_all_vars4, data=all_relev_vars)
summary(x4)

x5<-lm(time_trend_all_vars5, data=all_relev_vars)
summary(x5)

x6<-lm(time_trend_all_vars6, data=all_relev_vars)
summary(x6)

x7<-lm(time_trend_all_vars7, data=all_relev_vars)
summary(x7)

x8<-lm(time_trend_all_vars8, data=all_relev_vars)
summary(x8)

x9<-lm(time_trend_all_vars9, data=all_relev_vars)
summary(x9)

x10<-lm(time_trend_all_vars10, data=all_relev_vars)
summary(x10)

x11<-lm(time_trend_all_vars11, data=all_relev_vars)
summary(x11)

df_x<-broom::tidy(x)
df_x$term[grepl( "bri_week_" , df_x$term)]
df_x<-add_row(df_x,term="bri_week_10",estimate=0)
df_x<-subset(df_x, grepl( "bri_week_" , df_x$term))
df_x<-df_x[order(df_x$term),]

#Step12: Preparing GRAPH
x<-seq(-10, 42, by = 1)
length(x)
df_x$observation<-NA
df_x$observation <- x

df_x$observation<-as.character(df_x$observation)
selection<-df_x$observation[!stringr::str_detect(df_x$observation, '-')]
df_x$observation[!stringr::str_detect(df_x$observation, '-')]<-paste("+", selection, sep = "")
df_x$t<-paste("Networks x \n", "(t", df_x$observation, ")", sep = "")
df_x$t[df_x$t=="Networks x \n(t+0)"]<-"Network x t"

df_x$week<-NA
df_x$week<-readr::parse_number(df_x$term)
df_x$model<-"Bridging Networks"


###########################
#Zoom Event Study Bridging#
###########################

df_x_small<-subset(df_x, week>=(5) & week<=(16))

irislabs1<-seq(5,16,1)
irislabs2<-unique(df_x_small[as.numeric(df_x_small$week) %in% irislabs1,]$t)
length(irislabs2)

coef_zoom_x<-ggplot(df_x_small, aes(x=week, y=estimate))+
  geom_point(size = 2,
             position=position_dodge(width=0.4))+
  geom_errorbar(aes(ymin = estimate-1.96*std.error, ymax = estimate+1.96*std.error), 
                size = 1, width = 0,
                position=position_dodge(width=0.4))+
  scale_x_continuous(breaks = irislabs1,
                     labels = irislabs1,
                     sec.axis = sec_axis(~.,
                                         breaks =irislabs1,
                                         labels = irislabs2))+
  geom_hline(yintercept = 0) +
  geom_vline(xintercept=10,color="red")+
  #scale_x_continuous(name = "Bandwidth (km)", breaks=seq(5,85,5)) +
  labs(x = "Week", y = "Point Estimate")+
  theme_bw() +
  theme(axis.text.x = element_text(size=12),
        axis.text.x.top  = element_text(size=10, angle = 45, margin=margin(5,5,10,5,"pt")),
        axis.text.y = element_text(size=12),
        axis.title=element_text(size=12),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))
#plot_ex_coef2<-reposition_legend(plot_ex_coef, 'bottom left')
ggsave(coef_zoom_x, file = "paper/graphs/figure_a10k.jpg", 
       height = 15, width = 20, 
       units = "cm", dpi = 200)



###########################
#Controlling for pct_turnout#
###########################

df_x1<-broom::tidy(x1)
df_x1$term[grepl( "bri_week_" , df_x1$term)]
df_x1<-add_row(df_x1,term="bri_week_10",estimate=0)
df_x1<-subset(df_x1, grepl( "bri_week_" , df_x1$term))
df_x1<-df_x1[order(df_x1$term),]

#Step12: Preparing GRAPH
x<-seq(-10, 42, by = 1)
df_x1$observation<-NA
df_x1$observation <- x

df_x1$observation<-as.character(df_x1$observation)
selection<-df_x1$observation[!stringr::str_detect(df_x1$observation, '-')]
df_x1$observation[!stringr::str_detect(df_x1$observation, '-')]<-paste("+", selection, sep = "")
df_x1$t<-paste("Networks x \n", "(t", df_x1$observation, ")", sep = "")
df_x1$t[df_x1$t=="Networks x \n(t+0)"]<-"Network x t"

df_x1$week<-NA
df_x1$week<-readr::parse_number(df_x1$term)
df_x1$model<-"Bridging Networks"


###########################
#Zoom Event Study Bridging#
###########################

df_x1_small<-subset(df_x1, week>=(5) & week<=(16))

irislabs1<-seq(5,16,1)
irislabs2<-unique(df_x1_small[as.numeric(df_x1_small$week) %in% irislabs1,]$t)
length(irislabs2)

coef_zoom_x1<-ggplot(df_x1_small, aes(x=week, y=estimate))+
  geom_point(size = 2,
             position=position_dodge(width=0.4))+
  geom_errorbar(aes(ymin = estimate-1.96*std.error, ymax = estimate+1.96*std.error), 
                size = 1, width = 0,
                position=position_dodge(width=0.4))+
  scale_x_continuous(breaks = irislabs1,
                     labels = irislabs1,
                     sec.axis = sec_axis(~.,
                                         breaks =irislabs1,
                                         labels = irislabs2))+
  geom_hline(yintercept = 0) +
  geom_vline(xintercept=10,color="red")+
  #scale_x_continuous(name = "Bandwidth (km)", breaks=seq(5,85,5)) +
  labs(x = "Week", y = "Point Estimate")+
  theme_bw() +
  theme(axis.text.x = element_text(size=12),
        axis.text.x.top  = element_text(size=10, angle = 45, margin=margin(5,5,10,5,"pt")),
        axis.text.y = element_text(size=12),
        axis.title=element_text(size=12),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))
#plot_ex_coef2<-reposition_legend(plot_ex_coef, 'bottom left')
ggsave(coef_zoom_x1, file = "paper/graphs/figure_a10a.jpg", 
       height = 15, width = 20, 
       units = "cm", dpi = 200)




###########################
#Controlling for pop_total#
###########################

df_x2<-broom::tidy(x2)
df_x2$term[grepl( "bri_week_" , df_x2$term)]
df_x2<-add_row(df_x2,term="bri_week_10",estimate=0)
df_x2<-subset(df_x2, grepl( "bri_week_" , df_x2$term))
df_x2<-df_x2[order(df_x2$term),]

#Step12: Preparing GRAPH
x<-seq(-10, 42, by = 1)
length(x)
df_x2$observation<-NA
df_x2$observation <- x

df_x2$observation<-as.character(df_x2$observation)
selection<-df_x2$observation[!stringr::str_detect(df_x2$observation, '-')]
df_x2$observation[!stringr::str_detect(df_x2$observation, '-')]<-paste("+", selection, sep = "")
df_x2$t<-paste("Networks x \n", "(t", df_x2$observation, ")", sep = "")
df_x2$t[df_x2$t=="Networks x \n(t+0)"]<-"Network x t"

df_x2$week<-NA
df_x2$week<-readr::parse_number(df_x2$term)
df_x2$model<-"Bridging Networks"


###########################
#Zoom Event Study Bridging#
###########################

df_x2_small<-subset(df_x2, week>=(5) & week<=(16))

irislabs1<-seq(5,16,1)
irislabs2<-unique(df_x2_small[as.numeric(df_x2_small$week) %in% irislabs1,]$t)
length(irislabs2)

coef_zoom_x2<-ggplot(df_x2_small, aes(x=week, y=estimate))+
  geom_point(size = 2,
             position=position_dodge(width=0.4))+
  geom_errorbar(aes(ymin = estimate-1.96*std.error, ymax = estimate+1.96*std.error), 
                size = 1, width = 0,
                position=position_dodge(width=0.4))+
  scale_x_continuous(breaks = irislabs1,
                     labels = irislabs1,
                     sec.axis = sec_axis(~.,
                                         breaks =irislabs1,
                                         labels = irislabs2))+
  geom_hline(yintercept = 0) +
  geom_vline(xintercept=10,color="red")+
  #scale_x_continuous(name = "Bandwidth (km)", breaks=seq(5,85,5)) +
  labs(x = "Week", y = "Point Estimate")+
  theme_bw() +
  theme(axis.text.x = element_text(size=12),
        axis.text.x.top  = element_text(size=10, angle = 45, margin=margin(5,5,10,5,"pt")),
        axis.text.y = element_text(size=12),
        axis.title=element_text(size=12),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))
#plot_ex_coef2<-reposition_legend(plot_ex_coef, 'bottom left')
ggsave(coef_zoom_x2, file = "paper/graphs/figure_a10b.jpg", 
       height = 15, width = 20, 
       units = "cm", dpi = 200)


###########################
#Controlling for log_gdp#
###########################

df_x3<-broom::tidy(x3)
df_x3$term[grepl( "bri_week_" , df_x3$term)]
df_x3<-add_row(df_x3,term="bri_week_10",estimate=0)
df_x3<-subset(df_x3, grepl( "bri_week_" , df_x3$term))
df_x3<-df_x3[order(df_x3$term),]

#Step12: Preparing GRAPH
x<-seq(-10, 42, by = 1)
df_x3$observation<-NA
df_x3$observation <- x

df_x3$observation<-as.character(df_x3$observation)
selection<-df_x3$observation[!stringr::str_detect(df_x3$observation, '-')]
df_x3$observation[!stringr::str_detect(df_x3$observation, '-')]<-paste("+", selection, sep = "")
df_x3$t<-paste("Networks x \n", "(t", df_x3$observation, ")", sep = "")
df_x3$t[df_x3$t=="Networks x \n(t+0)"]<-"Network x t"

df_x3$week<-NA
df_x3$week<-readr::parse_number(df_x3$term)
df_x3$model<-"Bridging Networks"


###########################
#Zoom Event Study Bridging#
###########################

df_x3_small<-subset(df_x3, week>=(5) & week<=(16))

irislabs1<-seq(5,16,1)
irislabs2<-unique(df_x3_small[as.numeric(df_x3_small$week) %in% irislabs1,]$t)
length(irislabs2)

coef_zoom_x3<-ggplot(df_x3_small, aes(x=week, y=estimate))+
  geom_point(size = 2,
             position=position_dodge(width=0.4))+
  geom_errorbar(aes(ymin = estimate-1.96*std.error, ymax = estimate+1.96*std.error), 
                size = 1, width = 0,
                position=position_dodge(width=0.4))+
  scale_x_continuous(breaks = irislabs1,
                     labels = irislabs1,
                     sec.axis = sec_axis(~.,
                                         breaks =irislabs1,
                                         labels = irislabs2))+
  geom_hline(yintercept = 0) +
  geom_vline(xintercept=10,color="red")+
  #scale_x_continuous(name = "Bandwidth (km)", breaks=seq(5,85,5)) +
  labs(x = "Week", y = "Point Estimate")+
  theme_bw() +
  theme(axis.text.x = element_text(size=12),
        axis.text.x.top  = element_text(size=10, angle = 45, margin=margin(5,5,10,5,"pt")),
        axis.text.y = element_text(size=12),
        axis.title=element_text(size=12),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))
#plot_ex_coef2<-reposition_legend(plot_ex_coef, 'bottom left')
ggsave(coef_zoom_x3, file = "paper/graphs/figure_a10c.jpg", 
       height = 15, width = 20, 
       units = "cm", dpi = 200)


#########################
#Controlling for pop_den#
#########################

df_x4<-broom::tidy(x4)
df_x4$term[grepl( "bri_week_" , df_x4$term)]
df_x4<-add_row(df_x4,term="bri_week_10",estimate=0)
df_x4<-subset(df_x4, grepl( "bri_week_" , df_x4$term))
df_x4<-df_x4[order(df_x4$term),]

#Step12: Preparing GRAPH
x<-seq(-10, 42, by = 1)
df_x4$observation<-NA
df_x4$observation <- x

df_x4$observation<-as.character(df_x4$observation)
selection<-df_x4$observation[!stringr::str_detect(df_x4$observation, '-')]
df_x4$observation[!stringr::str_detect(df_x4$observation, '-')]<-paste("+", selection, sep = "")
df_x4$t<-paste("Networks x \n", "(t", df_x4$observation, ")", sep = "")
df_x4$t[df_x4$t=="Networks x \n(t+0)"]<-"Network x t"

df_x4$week<-NA
df_x4$week<-readr::parse_number(df_x4$term)
df_x4$model<-"Bridging Networks"


###########################
#Zoom Event Study Bridging#
###########################

df_x4_small<-subset(df_x4, week>=(5) & week<=(16))

irislabs1<-seq(5,16,1)
irislabs2<-unique(df_x4_small[as.numeric(df_x4_small$week) %in% irislabs1,]$t)
length(irislabs2)

coef_zoom_x4<-ggplot(df_x4_small, aes(x=week, y=estimate))+
  geom_point(size = 2,
             position=position_dodge(width=0.4))+
  geom_errorbar(aes(ymin = estimate-1.96*std.error, ymax = estimate+1.96*std.error), 
                size = 1, width = 0,
                position=position_dodge(width=0.4))+
  scale_x_continuous(breaks = irislabs1,
                     labels = irislabs1,
                     sec.axis = sec_axis(~.,
                                         breaks =irislabs1,
                                         labels = irislabs2))+
  geom_hline(yintercept = 0) +
  geom_vline(xintercept=10,color="red")+
  #scale_x_continuous(name = "Bandwidth (km)", breaks=seq(5,85,5)) +
  labs(x = "Week", y = "Point Estimate")+
  theme_bw() +
  theme(axis.text.x = element_text(size=12),
        axis.text.x.top  = element_text(size=10, angle = 45, margin=margin(5,5,10,5,"pt")),
        axis.text.y = element_text(size=12),
        axis.title=element_text(size=12),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))
#plot_ex_coef2<-reposition_legend(plot_ex_coef, 'bottom left')
ggsave(coef_zoom_x4, file = "paper/graphs/figure_a10d.jpg", 
       height = 15, width = 20, 
       units = "cm", dpi = 200)


#############################
#Controlling for pct_college#
#############################

df_x5<-broom::tidy(x5)
df_x5$term[grepl( "bri_week_" , df_x5$term)]
df_x5<-add_row(df_x5,term="bri_week_10",estimate=0)
df_x5<-subset(df_x5, grepl( "bri_week_" , df_x5$term))
df_x5<-df_x5[order(df_x5$term),]

#Step12: Preparing GRAPH
x<-seq(-10, 42, by = 1)
df_x5$observation<-NA
df_x5$observation <- x

df_x5$observation<-as.character(df_x5$observation)
selection<-df_x5$observation[!stringr::str_detect(df_x5$observation, '-')]
df_x5$observation[!stringr::str_detect(df_x5$observation, '-')]<-paste("+", selection, sep = "")
df_x5$t<-paste("Networks x \n", "(t", df_x5$observation, ")", sep = "")
df_x5$t[df_x5$t=="Networks x \n(t+0)"]<-"Network x t"

df_x5$week<-NA
df_x5$week<-readr::parse_number(df_x5$term)
df_x5$model<-"Bridging Networks"


###########################
#Zoom Event Study Bridging#
###########################

df_x5_small<-subset(df_x5, week>=(5) & week<=(16))

irislabs1<-seq(5,16,1)
irislabs2<-unique(df_x5_small[as.numeric(df_x5_small$week) %in% irislabs1,]$t)

coef_zoom_x5<-ggplot(df_x5_small, aes(x=week, y=estimate))+
  geom_point(size = 2,
             position=position_dodge(width=0.4))+
  geom_errorbar(aes(ymin = estimate-1.96*std.error, ymax = estimate+1.96*std.error), 
                size = 1, width = 0,
                position=position_dodge(width=0.4))+
  scale_x_continuous(breaks = irislabs1,
                     labels = irislabs1,
                     sec.axis = sec_axis(~.,
                                         breaks =irislabs1,
                                         labels = irislabs2))+
  geom_hline(yintercept = 0) +
  geom_vline(xintercept=10,color="red")+
  #scale_x_continuous(name = "Bandwidth (km)", breaks=seq(5,85,5)) +
  labs(x = "Week", y = "Point Estimate")+
  theme_bw() +
  theme(axis.text.x = element_text(size=12),
        axis.text.x.top  = element_text(size=10, angle = 45, margin=margin(5,5,10,5,"pt")),
        axis.text.y = element_text(size=12),
        axis.title=element_text(size=12),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))
#plot_ex_coef2<-reposition_legend(plot_ex_coef, 'bottom left')
ggsave(coef_zoom_x5, file = "paper/graphs/figure_a10e.jpg", 
       height = 15, width = 20, 
       units = "cm", dpi = 200)



##############################
#Controlling for pct_above_60#
##############################

df_x6<-broom::tidy(x6)
df_x6$term[grepl( "bri_week_" , df_x6$term)]
df_x6<-add_row(df_x6,term="bri_week_10",estimate=0)
df_x6<-subset(df_x6, grepl( "bri_week_" , df_x6$term))
df_x6<-df_x6[order(df_x6$term),]

#Step12: Preparing GRAPH
x<-seq(-10, 42, by = 1)
df_x6$observation<-NA
df_x6$observation <- x

df_x6$observation<-as.character(df_x6$observation)
selection<-df_x6$observation[!stringr::str_detect(df_x6$observation, '-')]
df_x6$observation[!stringr::str_detect(df_x6$observation, '-')]<-paste("+", selection, sep = "")
df_x6$t<-paste("Networks x \n", "(t", df_x6$observation, ")", sep = "")
df_x6$t[df_x6$t=="Networks x \n(t+0)"]<-"Network x t"

df_x6$week<-NA
df_x6$week<-readr::parse_number(df_x6$term)
df_x6$model<-"Bridging Networks"


###########################
#Zoom Event Study Bridging#
###########################

df_x6_small<-subset(df_x6, week>=(5) & week<=(16))

irislabs1<-seq(1,53,1)
irislabs2<-unique(df_x6[as.numeric(df_x6$week) %in% irislabs1,]$t)

coef_zoom_x6<-ggplot(df_x6_small, aes(x=week, y=estimate))+
  geom_point(size = 2,
             position=position_dodge(width=0.4))+
  geom_errorbar(aes(ymin = estimate-1.96*std.error, ymax = estimate+1.96*std.error), 
                size = 1, width = 0,
                position=position_dodge(width=0.4))+
  scale_x_continuous(breaks = irislabs1,
                     labels = irislabs1,
                     sec.axis = sec_axis(~.,
                                         breaks =irislabs1,
                                         labels = irislabs2))+
  geom_hline(yintercept = 0) +
  geom_vline(xintercept=10,color="red")+
  #scale_x_continuous(name = "Bandwidth (km)", breaks=seq(5,85,5)) +
  labs(x = "Week", y = "Point Estimate")+
  theme_bw() +
  theme(axis.text.x = element_text(size=12),
        axis.text.x.top  = element_text(size=10, angle = 45, margin=margin(5,5,10,5,"pt")),
        axis.text.y = element_text(size=12),
        axis.title=element_text(size=12),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))
#plot_ex_coef2<-reposition_legend(plot_ex_coef, 'bottom left')
ggsave(coef_zoom_x6, file = "paper/graphs/figure_a10f.jpg", 
       height = 15, width = 20, 
       units = "cm", dpi = 200)



###########################
#Controlling for pct_sec#
###########################

df_x7<-broom::tidy(x7)
df_x7$term[grepl( "bri_week_" , df_x7$term)]
df_x7<-add_row(df_x7,term="bri_week_10",estimate=0)
df_x7<-subset(df_x7, grepl( "bri_week_" , df_x7$term))
df_x7<-df_x7[order(df_x7$term),]

#Step12: Preparing GRAPH
x<-seq(-10, 42, by = 1)
df_x7$observation<-NA
df_x7$observation <- x

df_x7$observation<-as.character(df_x7$observation)
selection<-df_x7$observation[!stringr::str_detect(df_x7$observation, '-')]
df_x7$observation[!stringr::str_detect(df_x7$observation, '-')]<-paste("+", selection, sep = "")
df_x7$t<-paste("Networks x \n", "(t", df_x7$observation, ")", sep = "")
df_x7$t[df_x7$t=="Networks x \n(t+0)"]<-"Network x t"

df_x7$week<-NA
df_x7$week<-readr::parse_number(df_x7$term)
df_x7$model<-"Bridging Networks"


###########################
#Zoom Event Study Bridging#
###########################

df_x7_small<-subset(df_x7, week>=(5) & week<=(16))

irislabs1<-seq(5,16,1)
irislabs2<-unique(df_x7_small[as.numeric(df_x7_small$week) %in% irislabs1,]$t)
length(irislabs2)

coef_zoom_x7<-ggplot(df_x7_small, aes(x=week, y=estimate))+
  geom_point(size = 2,
             position=position_dodge(width=0.4))+
  geom_errorbar(aes(ymin = estimate-1.96*std.error, ymax = estimate+1.96*std.error), 
                size = 1, width = 0,
                position=position_dodge(width=0.4))+
  scale_x_continuous(breaks = irislabs1,
                     labels = irislabs1,
                     sec.axis = sec_axis(~.,
                                         breaks =irislabs1,
                                         labels = irislabs2))+
  geom_hline(yintercept = 0) +
  geom_vline(xintercept=10,color="red")+
  #scale_x_continuous(name = "Bandwidth (km)", breaks=seq(5,85,5)) +
  labs(x = "Week", y = "Point Estimate")+
  theme_bw() +
  theme(axis.text.x = element_text(size=12),
        axis.text.x.top  = element_text(size=10, angle = 45, margin=margin(5,5,10,5,"pt")),
        axis.text.y = element_text(size=12),
        axis.title=element_text(size=12),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))
#plot_ex_coef2<-reposition_legend(plot_ex_coef, 'bottom left')
ggsave(coef_zoom_x7, file = "paper/graphs/figure_a10g.jpg", 
       height = 15, width = 20, 
       units = "cm", dpi = 200)



###########################
#Controlling for pct_manuf#
###########################

df_x8<-broom::tidy(x8)
df_x8$term[grepl( "bri_week_" , df_x8$term)]
df_x8<-add_row(df_x8,term="bri_week_10",estimate=0)
df_x8<-subset(df_x8, grepl( "bri_week_" , df_x8$term))
df_x8<-df_x8[order(df_x8$term),]

#Step12: Preparing GRAPH
x<-seq(-10, 42, by = 1)
df_x8$observation<-NA
df_x8$observation <- x

df_x8$observation<-as.character(df_x8$observation)
selection<-df_x8$observation[!stringr::str_detect(df_x8$observation, '-')]
df_x8$observation[!stringr::str_detect(df_x8$observation, '-')]<-paste("+", selection, sep = "")
df_x8$t<-paste("Networks x \n", "(t", df_x8$observation, ")", sep = "")
df_x8$t[df_x8$t=="Networks x \n(t+0)"]<-"Network x t"

df_x8$week<-NA
df_x8$week<-readr::parse_number(df_x8$term)
df_x8$model<-"Bridging Networks"


###########################
#Zoom Event Study Bridging#
###########################

df_x8_small<-subset(df_x8, week>=(5) & week<=(16))

irislabs1<-seq(1,53,1)
irislabs2<-unique(df_x8[as.numeric(df_x8$week) %in% irislabs1,]$t)

coef_zoom_x8<-ggplot(df_x8_small, aes(x=week, y=estimate))+
  geom_point(size = 2,
             position=position_dodge(width=0.4))+
  geom_errorbar(aes(ymin = estimate-1.96*std.error, ymax = estimate+1.96*std.error), 
                size = 1, width = 0,
                position=position_dodge(width=0.4))+
  scale_x_continuous(breaks = irislabs1,
                     labels = irislabs1,
                     sec.axis = sec_axis(~.,
                                         breaks =irislabs1,
                                         labels = irislabs2))+
  geom_hline(yintercept = 0) +
  geom_vline(xintercept=10,color="red")+
  #scale_x_continuous(name = "Bandwidth (km)", breaks=seq(5,85,5)) +
  labs(x = "Week", y = "Point Estimate")+
  theme_bw() +
  theme(axis.text.x = element_text(size=12),
        axis.text.x.top  = element_text(size=10, angle = 45, margin=margin(5,5,10,5,"pt")),
        axis.text.y = element_text(size=12),
        axis.title=element_text(size=12),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))
#plot_ex_coef2<-reposition_legend(plot_ex_coef, 'bottom left')
ggsave(coef_zoom_x8, file = "paper/graphs/figure_a10h.jpg", 
       height = 15, width = 20, 
       units = "cm", dpi = 200)


###########################
#Controlling for police_den#
###########################

df_x9<-broom::tidy(x9)
df_x9$term[grepl( "bri_week_" , df_x9$term)]
df_x9<-add_row(df_x9,term="bri_week_10",estimate=0)
df_x9<-subset(df_x9, grepl( "bri_week_" , df_x9$term))
df_x9<-df_x9[order(df_x9$term),]

#Step12: Preparing GRAPH
x<-seq(-10, 42, by = 1)
df_x9$observation<-NA
df_x9$observation <- x

df_x9$observation<-as.character(df_x9$observation)
selection<-df_x9$observation[!stringr::str_detect(df_x9$observation, '-')]
df_x9$observation[!stringr::str_detect(df_x9$observation, '-')]<-paste("+", selection, sep = "")
df_x9$t<-paste("Networks x \n", "(t", df_x9$observation, ")", sep = "")
df_x9$t[df_x9$t=="Networks x \n(t+0)"]<-"Network x t"

df_x9$week<-NA
df_x9$week<-readr::parse_number(df_x9$term)
df_x9$model<-"Bridging Networks"


###########################
#Zoom Event Study Bridging#
###########################

df_x9_small<-subset(df_x9, week>=(5) & week<=(16))

irislabs1<-seq(5,16,1)
irislabs2<-unique(df_x9_small[as.numeric(df_x9_small$week) %in% irislabs1,]$t)
length(irislabs2)

coef_zoom_x9<-ggplot(df_x9_small, aes(x=week, y=estimate))+
  geom_point(size = 2,
             position=position_dodge(width=0.4))+
  geom_errorbar(aes(ymin = estimate-1.96*std.error, ymax = estimate+1.96*std.error), 
                size = 1, width = 0,
                position=position_dodge(width=0.4))+
  scale_x_continuous(breaks = irislabs1,
                     labels = irislabs1,
                     sec.axis = sec_axis(~.,
                                         breaks =irislabs1,
                                         labels = irislabs2))+
  geom_hline(yintercept = 0) +
  geom_vline(xintercept=10,color="red")+
  #scale_x_continuous(name = "Bandwidth (km)", breaks=seq(5,85,5)) +
  labs(x = "Week", y = "Point Estimate")+
  theme_bw() +
  theme(axis.text.x = element_text(size=12),
        axis.text.x.top  = element_text(size=10, angle = 45, margin=margin(5,5,10,5,"pt")),
        axis.text.y = element_text(size=12),
        axis.title=element_text(size=12),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))
#plot_ex_coef2<-reposition_legend(plot_ex_coef, 'bottom left')
ggsave(coef_zoom_x9, file = "paper/graphs/figure_a10j.jpg", 
       height = 15, width = 20, 
       units = "cm", dpi = 200)



###########################
#Controlling for east_ger#
###########################

df_x10<-broom::tidy(x10)
df_x10$term[grepl( "east_ger_week_" , df_x10$term)]
df_x10<-add_row(df_x10,term="east_ger_week_10",estimate=0)
df_x10<-subset(df_x10, grepl( "east_ger_" , df_x10$term))
df_x10<-df_x10[order(df_x10$term),]

#Step12: Preparing GRAPH
x<-seq(-10, 42, by = 1)
df_x10$observation<-NA
df_x10$observation <- x

df_x10$observation<-as.character(df_x10$observation)
selection<-df_x10$observation[!stringr::str_detect(df_x10$observation, '-')]
df_x10$observation[!stringr::str_detect(df_x10$observation, '-')]<-paste("+", selection, sep = "")
df_x10$t<-paste("Networks x \n", "(t", df_x10$observation, ")", sep = "")
df_x10$t[df_x10$t=="Networks x \n(t+0)"]<-"Network x t"

df_x10$week<-NA
df_x10$week<-readr::parse_number(df_x10$term)
df_x10$model<-"Bridging Networks"

###########################
#Zoom Event Study Bridging#
###########################

df_x10_small<-subset(df_x10, week>=(5) & week<=(16))

irislabs1<-seq(5,16,1)
irislabs2<-unique(df_x10_small[as.numeric(df_x10_small$week) %in% irislabs1,]$t)

coef_zoom_x10<-ggplot(df_x10_small, aes(x=week, y=estimate))+
  geom_point(size = 2,
             position=position_dodge(width=0.4))+
  geom_errorbar(aes(ymin = estimate-1.96*std.error, ymax = estimate+1.96*std.error), 
                size = 1, width = 0,
                position=position_dodge(width=0.4))+
  scale_x_continuous(breaks = irislabs1,
                     labels = irislabs1,
                     sec.axis = sec_axis(~.,
                                         breaks =irislabs1,
                                         labels = irislabs2))+
  geom_hline(yintercept = 0) +
  geom_vline(xintercept=10,color="red")+
  #scale_x_continuous(name = "Bandwidth (km)", breaks=seq(5,85,5)) +
  labs(x = "Week", y = "Point Estimate")+
  theme_bw() +
  theme(axis.text.x = element_text(size=12),
        axis.text.x.top  = element_text(size=10, angle = 45, margin=margin(5,5,10,5,"pt")),
        axis.text.y = element_text(size=12),
        axis.title=element_text(size=12),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))
#plot_ex_coef2<-reposition_legend(plot_ex_coef, 'bottom left')
ggsave(coef_zoom_x10, file = "paper/graphs/figure_a10i.jpg", 
       height = 15, width = 20, 
       units = "cm", dpi = 200)



###########################
#Controlling for everything#
###########################

df_x11<-broom::tidy(x11)
df_x11$term[grepl( "bri_week_" , df_x11$term)]
df_x11<-add_row(df_x11,term="bri_week_10",estimate=0)
df_x11<-subset(df_x11, grepl( "bri_week_" , df_x11$term))
df_x11<-df_x11[order(df_x11$term),]

#Step12: Preparing GRAPH
x<-seq(-10, 42, by = 1)
df_x11$observation<-NA
df_x11$observation <- x

df_x11$observation<-as.character(df_x11$observation)
selection<-df_x11$observation[!stringr::str_detect(df_x11$observation, '-')]
df_x11$observation[!stringr::str_detect(df_x11$observation, '-')]<-paste("+", selection, sep = "")
df_x11$t<-paste("Networks x \n", "(t", df_x11$observation, ")", sep = "")
df_x11$t[df_x11$t=="Networks x \n(t+0)"]<-"Network x t"

df_x11$week<-NA
df_x11$week<-readr::parse_number(df_x11$term)
df_x11$model<-"Bridging Networks"

###########################
#Zoom Event Study Bridging#
###########################

df_x11_small<-subset(df_x11, week>=(5) & week<=(16))

irislabs1<-seq(5,16,1)
irislabs2<-unique(df_x11_small[as.numeric(df_x11_small$week) %in% irislabs1,]$t)
length(irislabs2)

coef_zoom_x11<-ggplot(df_x11_small, aes(x=week, y=estimate))+
  geom_point(size = 2,
             position=position_dodge(width=0.4))+
  geom_errorbar(aes(ymin = estimate-1.96*std.error, ymax = estimate+1.96*std.error), 
                size = 1, width = 0,
                position=position_dodge(width=0.4))+
  scale_x_continuous(breaks = irislabs1,
                     labels = irislabs1,
                     sec.axis = sec_axis(~.,
                                         breaks =irislabs1,
                                         labels = irislabs2))+
  geom_hline(yintercept = 0) +
  geom_vline(xintercept=10,color="red")+
  #scale_x_continuous(name = "Bandwidth (km)", breaks=seq(5,85,5)) +
  labs(x = "Week", y = "Point Estimate")+
  theme_bw() +
  theme(axis.text.x = element_text(size=12),
        axis.text.x.top  = element_text(size=10, angle = 45, margin=margin(5,5,10,5,"pt")),
        axis.text.y = element_text(size=12),
        axis.title=element_text(size=12),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))
#plot_ex_coef2<-reposition_legend(plot_ex_coef, 'bottom left')
ggsave(coef_zoom_x11, file = "paper/graphs/figure_a10k.jpg", 
       height = 15, width = 20, 
       units = "cm", dpi = 200)

